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Abstract 

Multispecies reaction-diffusion systems, for which the time evolution 
equation of correlation functions become a closed set, are considered. A 
formal solution for the average densities is found. Some special inter- 
actions and the exact time dependence of the average densities in these 
cases are also studied. For the general case, the large time behaviour of 
the average densities has also been obtained. 



1 Introduction 

In recent years, reaction-diffusion systems have been studied by many peo- 
ple. As mean field techniques, generally do not give correct results for low 
dimensional systems, people are motivated to study stochastic models in low 
dimensions. Moreover, solving one dimensional systems should in principle be 
easier. Exact results for some models in a one-dimensional lattice have been 
obtained, for example in [1-10]. 

Different methods have been used to study these models, including analyt- 
ical and asymptotic methods, mean field methods, and large-scale numerical 
methods. Systems with more than one species have also been studied [11-23]. 
Most of the arguments are based on simulation results. There are, however, 
some exact results as well ( |l8|, ||(], ^] for example). 

In p3| , a 10-parameter family of stochastic models has been studied. In 
these models, the fc-point equal time correlation functions (niUj ■ ■ ■ n^) satisfy 
linear differential equations involving no higher-order correlations. These linear 
equations for the average density (rii) has been solved. But these set of equa- 
tions can not be solved easily for higher order correlation functions. We have 
generalized the same idea to multi-species models. We have considered general 
reaction diffusion processes of multi-species in one dimension with two-site in- 
teraction. We have obtained the conditions the Hamiltonian should satisfy in 
order to give rise to closed set of time evolution equation for correlation func- 
tions. The set of equations for average densities can be written in terms of four 
matrices. The time evolution equation for more-point functions, besides these 
four matrices, generally depend explicitly on the elements of the Hamiltonian 
, and generally can not be solved easily. These matrices are not determined 
uniquely from the Hamiltonian: there is a kind of gauge transformation one 
can apply on them which of course, does not change the evolution equation. 
A formal solution for average densities of different species is found. For some 
special choices of the four matrices we also give the explicit form of interactions 
and the exact time dependence of average densities. At the end, we study the 
large time behaviour of the average densities of different species for the general 
case. 

2 A brief review of linear stochastic systems 

To fix the notation used in this article, here we briefly review the already well 
known formalism of linear stochastic systems. The master equation for P(a, t) 



where lo(t — > a) is the transition rate from the configuration r to a. Introducing 
the state vector 



is 




(1) 




(2) 



a 



l 



where the summation runs over all possible states of the system, one can write 
the above equation in the form 

±\P)=H\P), (3) 
where the matrix elements of H. are 

{u\H\t) = w(r ->■ a), r ^ a, 

(<j\H\a) = -J2^^r). (4) 

The basis {(oj} is dual to that is 

(<t|t) = S aT . (5) 

The operator Ti. is called a Hamiltonian, and it is not necessarily hermitian. 
Conservation of probability, 

£P((7,t)=l, (6) 
a 

shows that 

(S\H = 0, (7) 

where 

(S|=5>|. (8) 



So, the sum of each column of H., as a matrix, should be zero. As (S\ is a 
left eigenvector of 7i with zero eigenvalue, 7i has at least one right eigenvector 
with zero eigenvalue. This state corresponds to the steady state distribution of 
the system and it does not evolve in time. If the zero eigenvalue is degenerate, 
the steady state is not unique. The transition rates are non-negative, so the 
off-diagonal elements of the matrix 7i are non-negative. Therefore, if a matrix 
H has the following properties, 

(s\n = o, 

(a\H\r) > 0, W 

then it can be considered as the generator of a stochastic process. It can be 
proved that the real part of the eigenvalues of any matrix with the above con- 
ditions is less than or equal to zero. 

The dynamics of the state vectors (||) is given by 

\P(t))=exp(m)\P(fl)), (10) 

and the expectation value of an observable O is 

(0)(t) =Y,0{cT)P{a,t) = (S\Oexp(tH)\P(0)). (11) 
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3 Models leading to closed set of evolution equa- 
tions 



The models which we address are multispecies reaction-diffusion models. That 
is, each site is a vacancy or has one particle. There are several kinds of par- 
ticles, but at any time at most one kind can be present at each site. Suppose 
the interaction is between nearest neighbors, and the system is translationally 
invariant. 

L 

U = ^£f i)i+ i. (12) 

i=l 

The number of sites is L and the number of possible states in a site is TV; 
different states of each site are denoted by A a , a = 1, • • • N, where one of the 
states is vacancy. Introducing nf as the number operator of A a particle in the 
site i, we have 

£< = !■ (13) 

a=l 

The average number density of the particle A a in the site i at the time t is 

(nf) = (S\nf\P(t)) (14) 

where \P(t)) := cxp(tH) \P (0)) represents the state of the system at the time t, 

(S| = <s|®"-®<s|, (15) 

V v ' 

L 

and 

(s\ :=(1 1 ...1). (16) 

S v ' 

N 

So, the time evolution of (nf ) is given by 

±(nf) = (S\nfH\P(t)). (17) 

The only terms of the Hamiltonian TL which are relevant in the above equation 
are -ff^i+i and Hi-i^. The result of acting any matrix Q on the ket (s\ is 
equivalent to acting the diagonal matrix Q on the same ket, provided each 
diagonal element of the matrix Q is the sum of all elements of the corresponding 
column in the matrix Q. So, the action of (1 £g> n a )H and (n a <g> \)H on (s\ <g> (s\ 
are equivalent to the action of two diagonal matrices on (s\ ® (s\. We use the 
notation ~, for the equivalent action on (s\ ® (s\. 

[\®n a )H - Yj A Pi nfi ® nl 
Pi 
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(n a (g>l)H ~ nP ® nl - (18) 

where Ap and „4^j 7 are as the following 

A 
A 

Then, equation (|lj]) takes the following form 

«) =X)^<"f-i»7> + -^>fn? +1 >. (20) 

Hi 

Generally, in the time evolution equation of (n a ) the two-point functions (n^n 7 ) 
appear. Using ([l3]), one can see that iff A and A satisfy the following equations, 
then the right hand side of the ( |20| ) can be expressed in terms of only one-point 
functions. 

AO. , AO. _ AO. _ AO. _ n 

A a Pl +A a NN -A% 7 -A a m = o. (21) 

These equations give 2(N — l) 3 constraints on the Hamiltonian, so adding the 
condition of stochasticity of H, we have 2(N — l) 3 + N 2 relations between the 
elements of H. The constraints (Ell) mean 



A a /-ia -not 

JKp^ — ljp_ — -D 7 _ 

A% = —Bp + D". (22) 

So, ( |l7| ) takes the form 

[ - {Bf + 5|)(nf ) + C^nf^) + D|<nf +1 )] . (23) 



N 



(3=1 

In the simplest case, the one-species, each site is vacant or occupied by only 
one kind of particles. Then, the matrices B, C, B, and D are two-dimensional. 
Using (|l3|), The equation for (hj) is 

(n\) = {-B\- B\ + B\ + B\) (nj) + (C\ - d) (n\_ x ) 

+ {D\ - D\) (nl +1 ) + ( - Bl - B\ + C\ + D\) . (24) 

This is a linear difference equation, of the kind obtained [^3| , and its solution 
can be expressed in terms of modified Bessel functions. 

The time evolution equation for two-point functions also can be obtained. 

d 



dt 



«<> = E[-( S 7+^)(">f) + C 7>7-i^) + ^(n7 + i^) 
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-(^ + ^)<rz>J)+^«nJ_ 1 ) 

+D^nfn] +1 )), \i-j\>l, (25) 



<l nfnL.) = T[- Banjul,) - S»7 +1 ) 



+C>7_ 1 nf +1 > +D^{nfn] +2 )] + ]T Hf^njn^). (26) 

7 A 

For more-point functions, one can deduce similar results. In fact, it is easy 
to show that if the evolution equations of one-point functions are closed, the 
evolution equation of n-point functions contain only n- and less-point functions. 
However, generally these set of equations can not be solved easily. 

4 Equivalent Hamiltonians regarding one-point 
functions, and gauge transformations 

Knowing B, C, B, and D does not determine the Hamiltonian uniquely, but as 
it is seen from (^3|), the time evolution of one-point functions depends only on 
B, C, B, and D. The two- and more-point functions depend explicitly on the 
elements of H. So, different Hamiltonians may give same evolutions for (nf). 
Take two Hamiltonians H and H' . Defining 

R:=H-H', (27) 

E< = E< = ' ( 28 ) 

a 

these two Hamiltonians give rise to the same A and A. Regarding one-point 
functions (nf), these models are the same. So, we call these models, regarding 
one-point functions, equivalent. 

However, A and A do not determine B, C, B, and D uniquely. The stochas- 
tic condition 

5> 7 t = 0, (29) 

a/3 

results in some constraints on B, C, B, and D: 

Y^iCp - B°) = 

a. 

J2(-B% + D«) = 0, (30) 

OL 

So, the sum of all elements of any column of B (C) should be the same 

£ 0-? = $>£ = /. (31) 
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Then, the state (s| is the left eigenvector of B and C, with the same eigenvalue 
/. B and D have also the same property, of course with different eigenvalue g. 
Changing B and C according to the gauge transformation, 

C%^C ia p=C%-t or C' = C-\f)(s\ 

B$ - B'l = B$ - f a or B'=B-\f){s\ (32) 

does not change A . With a suitable choice of 

a 

the sum of the elements of any column of B or C can be set to zero. In this 
gauge, the eigenvalues of B and C for the eigenvector (s| will be zero. 



5 One-point functions 

To solve (HI), we introduce the vector Afk 

( H> \ 



A4 



\«)J 



Equation ( |23|) can then be written as 

A4 = -(B + B)M k + CAA fe _! + DN k+ i. 
Introducing the generating function G(z, t), 

oo 

G(z,i)=^A4(i)z fc , 

— oo 

one arrives at, 

G(z, t) = [ - (B + B) + z C + z- 1 D) G(z, t), 
the solution to which is 

G(z, t) = exp (t[-{B + B) + zC + z- 1 D])G(z, 0). 
A4(0' s are the coefficients of the Laurent expansion of G(z, t), so 



(34) 



(35) 

(36) 

(37) 
(38) 



A4(i) = ^- jr I dz z m - k - 1 exp(t[-(B + B) + zC + z- 1 D])Af m {0). (39) 



m— — 00 



This is the formal solution of the problem, which is of the form 

A4(*) = ^r fem (*yv m (o). (40) 



G 



5.1 Some special cases 

We now consider special choices for B, C, B, and D. 



B = 


\u)(b\, 


C = 


\u)(c\, 


B = 


\u)(b\, 


D = 


\u)(d\, 


\u) := ( 







5.1.1 The matrices B, C, B, and D are two dimensional (the single 
species case) 

We can use the gauge transformation to make (s\ the simultaneous null left- 
eigenvector of B, C, B, and D. In this gauge, one has 



(41) 

where 

(42) 

This means that it is orthogonal to (s\, and is a simultaneous right-eigenvector 
of B, C, B, and D. Using (|4l|), one can easily calculate the exponential in (|39|): 

(,tg(z) _ i 

exp (t[-(B + B) + z C + z- 1 D] =1+ \u){g{z)\, (43) 

where 

{g{z)\:=-{b\-{b\+z{c\+z-^dl (44) 

and 

g(z) := (g(z)\u). (45) 

Now take (t;| and \w) to be the left-eigenvector oi—B — B + C + D dual to \u) , 
and the right-eigenvector of — B — B + C + D dual to (s\, respectively. One can 
normalize these, so that 

(v\u) = 1, 

(s\w) = 1. (46) 
Of course, (v\ is orthogonal to \w). Then, 

_ P *9(z) _ 1 

exp (t[-(B+B)+z C+z' 1 D]) = e tg ^\u){v\ + \w){s\ + {g(z)\w) -— |u)(s|. 

g{z) 

(47) 

Acting this on Af m (0), and noting that 

(s\N m (0) - 1, (48) 

it is seen that 

oo 

N k (t) = \w)( 3 \N k (0) + — V <b dz z^-WWMivltfniQ) 

2m ^— ' J 



m— — OO 
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oo 

+ — V idz z^-'e^^iv^iO), (49) 
Zm *■ — ' / 

rn— — oo 



or 



1 /" 

(«|iV fc (t) = — £ i dzz^-^'^luX^IMnW. (50) 
This is equivalent to 



m— — oo 



5.1.2 C=pB, D = qB 
Using © 

(l-p)(*|B = 

means that p = 1 or (s|S = 0, and g = 
eigenvector of B and _B, then p = q = 
Now, using the definition of A 



■ (l-q)(s\B = 0. (51) 

= 1 or (s\B = 0. If (s\ is not the left null 
; 1. So, we will have B = C, and D = B. 



X 

A« =C«-C% = Y,H%. (52) 

A 

For a ^ /3, and all the terms in the right hand side summations in the 

above equations are reaction rates and should be non-negative, but the sum of 
the left hand sides is zero. So, 

C a = C a = f a i for 1 ^ a ^p. (53) 

All the elements of each row except the diagonal elements of C ( or B ) are the 
same. That is, 

C=\f)(s\ + C, (54) 

where C is some diagonal matrix. The fact that \s) is a left-eigenvector of C, 
shows that it should be a left-eigenvector of C as well. And this demands C 
to be proportional to the unit matrix. One can do the same arguments for B 
and D. So, after gauge transformation, 

C = B = ul, D = B = vl. (55) 

Although, the time evolution of average densities can be written in terms of 
B, C, B, and D , the Hamiltonian H is not uniquely be determined by these 
matrices. There exist different Hamiltonians which are equivalent, regarding 
one-point functions. 

x 

Y,H$ = A% = v{5«-8%). (56) 
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All the elements of the (3(3 column of H are zero. For a ^ (3, the elements of H 
satisfy 

E H a0 + H a0 + H % = U 

EO-/3A rrPa rr/3/3 _ 
H af3 + H al3 + H a0 - V 

E H^ + H^+Hf^-v. (57) 

In general, these sets of equations have several solutions, but for the one-species 
case, the reaction rates are the following 



A0 -> { AA u- A 12 (58) 



0A 


Al2 




AA 


u — 


A 


00 


u — 


A; 


A0 


A 2 i 




AA 


u — 


A : 


00 


u — 


A. 



®A—t\ AA v- A 21 (59) 

^21 

The above system, with no diffusion, has been studied in [ p4| . There, the n- 
point functions have been investigated. This solution can be generalized to the 
multispecies case. For a ^ (3 

A a Ap —> ApA a , A aP a,/3 = l---N 

A a Ap > A a A a , u — A aj g 

A a A fj -» « - A aP . (60) 

The only constraint is the non-negativeness of the reaction rates: 

u > A Q/3 > v> A Q(3 > 0. (61) 

This model has N(N— 1) + 2 free parameters. However, only the two parameters 
u and v appear in the time evolution equation of average densities: 

(hf) = -(« + «)«) + u{nU) + ««+!>• (62) 

As it is seen, dynamics of average densities of different particles decouple, and 
despite the complex interactions of the model, (nf )'s can be easily calculated. 
But in the time evolution of two-point functions A a/ g's appear as well. So, al- 
though models with different exchanging rates (A a/ g) and same initial conditions 
have the same average densities, their two-point functions generally are not the 
same. 
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5.1.3 B, B, C, D commute 

Generally, the gauge transformation do not preserve the commutation relation 
of B and C ( and that of B and D ). But if B and C commute, there is a gauge 
transformation which leaves the transformed B and C commuting. If we choose 
|/) to be a right eigenvector of B and C dual to (s|, that is 

B\f) = C\f)=f\f), (63) 

then B' := B — \f)(s\ and C := C - \f)(s\ commute. If 

(s\f) = f, (64) 

then (s| times B' and C will be zero. So, if B, C, B, and D commute with each 
other, there exists a suitable gauge transformation that makes their eigenvalue 
corresponding to (s| zero, while they remain commuting: 

(s\B = (s\C = (s\B = (s\D = 0, (65) 

Denote, the matrix which simultaneously diagonalize these four matrices by U , 
diagonalized matrices by primes, and their eigenvalues by b a , c a , b a , and d a , 
respectively. We have 

(n\B' = (0|C = (n\B' = (n\D' = o 

(fij = (s\U. (66) 
We take b N = c N = b N = d N = 0, and normalize (fi| and U so that 

= (0 ••■ 0), (67) 

and 

^2U a p = S Nfi . (68) 

a 

U will also diagonalize the exponential in (^). So we have 
1 00 r 

^ = 2^ ^ fdzz^-^expiti-B'-B' + zC' + z-'D'})^), (69) 



771— — OO 



where 

K(t) := U- l Af k (t). (70) 

The matrix in the argument of the exponential in (^9|) is diagonal, so the integral 
can be easily calculated: 

I := — I dz z™-^ 1 exp (t[-b a -b a + z c a + z' 1 d a }). (71) 
2m J 

Introducing w := sj^z, one arrives at 



I— (_ ) i t Jl I dw w m - k - 1 exp ( Vc a d a t(w + w- 1 )) , (72) 

c a 2ni 1 
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which can be written in terms of modified Bessel functions 
Then, 



AT fc (t)= X! f/diag|(^)^ e - t ^+^)l^„ i (2v / ^t)|L/- 1 AA m (0). (74) 



Note that the right-hand side of ( |73[ ) is 5k, m for a = N, since the iV-th eigen- 
value of B, C, B, and D is zero. 

One can start with four special diagonal matrices, and then construct the 
Hamiltonians with different reaction-diffusion rates. Not all diagonal matrices 
lead to physical stochastic models: negative reaction rates may be obtained. 
Considering the large time behaviour of average number densities, one can show 
that 

|Re(\/c«J«)| < Re(b a + b a ), (75) 

which also shows that 

Re(b a +b a )>0. (76) 
Now, we consider a special choice for U : 

U$ = S%-(l-5%)5%. (77) 

Then 





= 6/3 


-5%) 




= cp(#f) 


-Sn) 




= hi&t 


-8%) 




= d p (S$ 





Now, consider 



(78) 



A %i = E H £ = - ba5 i + c ^ + + (79) 

A 

For a 7^ 7 and a ^ /?, 

E^7>° E^>°- (80) 

A A 

So, taking (3,j^N and a = N, 

b 1 > cP. (81) 
The same reasoning is true for b 1 and dP : 

P > d p . (82) 
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Here too, similar to the previous example, the above choices for B, C, B, and 
D do not determine H uniquely. One particular solution for the reaction rates 
is 



For a ^ N 





' A a A N , 


Aatq, 




A N A a -► | 




d Q - 


Aat q 






6 Q - 


Aat q 






A a N 




A a A N -> | 






A a AT 






6 a - 


AqAT 



(83) 
(84) 



and, for a,f3^=N 



A a A N , bp — Ca — A a p 

/9j b a — dp — A a f3 (85) 

For a 7^ /3, the following reactions may also occur. For a < (3 



and for a > (3 



The constraint of non-negativeness of the reaction rates leads to 

c a < dp < c 7 a < /3 < 7 

< A a p < bp - Ca 

A a p <b a — dp 

< AjVa < da 

< A Na < b_ a 

< AaN < b a 

0<A aN < C a (88) 

5.1.4 type— change invariance 

Suppose B, C, B, and D have the property 

B a ptl= B l (89) 

and the same for the other three matrices. Note that the indices of these matrices 
are defined periodically, so that N + a, as an index, is equivalent to a. This 
is in fact a special case of commuting matrices discussed earlier. One can use 
([74|). To do so, one should know the simultaneous eigenvectors of B, C, B, 
and D, and their corresponding eigenvalues. It is not difficult to see that the 
eigenvectors are 

^-tWW < 90 » 



12 



The corresponding eigenvectors of B, for example, are 



(91) 



Finally, the matrix elements of the inverse of U are 



(IT 1 )? = -JLexp 



N 



i2ira{3\ 



N 



(92) 



These can be put directly in (74) 



6 Large time behaviour of average densities 

The large-time behaviour of the system is deduced through a steepest-descent 
analysis of the formal solution ([39]). One should consider the eigenvalues and 
the eigenvectors of the z-dependent matrix 

M(z) := —{B + B) + zC + z~ 1 D. (93) 

Denote the eigenvalues of this matrix by A Q (z). As for any value of z, the 
matrix M has {s\ as its left eigenvector corresponding to the eigenvalue zero, M 
will have a right eigenvector \w) dual to (s\. \w) is z-dependent, but one can 
normalize it so that 

(s\w(z)) = 1. (94) 

The fact that Af should not blow up at t ~ * oo assures that the real-part of the 
eigenvalues of M(z) are non-positive (at least for |z| = 1). If all of the other 
eigenvalues have negative real-parts, then at t — > oo only \w) survives. That is, 

ATfc(oo) = V /dzz m - fe -Xz))( S |AA m (0) 

m J 

= h(l)>- (95) 

We have used (s\Af m (0) = 1. This could also be obtained directly, using the evo- 
lution equation (|2^) , by setting equal to zero and assuming A4 independent 
of k. So, the final state of the system is the eigenvector of — (B + B) + C + D, 
corresponding to the eigenvalue zero. 

To investigate the next-to-leading term at t — * oo, consider the other eigen- 
values of M(z). Suppose that at z — Zq, X a is stationary. There may be more 
than one point having this property. So, we will have a set consisting of ZQ a 's. 
Each of these points corresponds to a stationary eigenvalue Aq . We choose that 
Zg a the corresponding eigenvalue of which has the largest real-part. Denote this 
point by zq, its corresponding stationary eigenvalue by Ao, and its corresponding 
right eigenvector by |i>o). The next -to-leading term in Af is then 

A/f ~z -V A °K>. (96) 
Note that zq is not necessarily a phase, its modulus may be different from 1. 
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